%%%%%%%%%% Function FindFermi %%%%%%%%%%
    function [Sx, Ef] = FindFermiOld(Energy, DOS1, DOS2, x1, Vx1, BackGround1, KT)
    % we use a coarser spatial grid here
    skip=20;
    x=x1(1:skip:end);
    Vx=Vx1(1:skip:end);
    BackGround=BackGround1(1:skip:end);
    Ne = length(Energy); 
    he = Energy(2)-Energy(1);
    Nx = length(Vx);
    nc = 1+Nx/2;                        % Even Grid Points
    hx = x(2)-x(1);
    S = zeros(1,Nx); 
    Ef = Energy(round((Ne)/2));     % Initial guess of fermi level (middle)
    Efmin = Energy(1);
    Efmax = Energy(end); 

    tic;
    for i = 1:30
        % Adding up charge at each grid point
        F = FermiFunction(Energy, Ef+Vx, KT); % 0.1 seconds
        S(1:nc-1) = DOS1*F(:,1:nc-1)*he; % 0.005 seconds
        S(nc:end) = DOS2*F(:,nc:end)*he; % 0.005 mseconds
        TotalCharge = sum(BackGround - S)*hx; %  0.000038 seconds
        % Binary Searching
        if TotalCharge > 0
            Efmin = Ef;
            Ef = (Ef + Efmax)/2;
        else
            Efmax = Ef;
            Ef = (Efmin + Ef)/2;
        end
    end;
    toc;
    F=FermiFunction(Energy, Ef+Vx1, KT);
    Nx1 = length(Vx1);
    nc1 = 1+(Nx1)/2;
    S1 = zeros(1,Nx1);
    S1(1:nc1-1) = DOS1*F(:,1:nc1-1)*he;
    S1(nc1:end) = DOS2*F(:,nc1:end)*he;
    Sx = BackGround1 - S1;
    end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%